"""
    BackgroundCheck

Background departure quality control comparing observations to model background
(first guess). This is the primary QC method in variational data assimilation.

# Mathematical Foundation

## Innovation
```
d = y - H(x_b)
```
where:
- y: observation value
- H: observation operator
- x_b: background state

## Total Error Variance
```
σ²_total = σ²_obs + σ²_bg = R + H*B*H^T
```
where:
- R: observation error covariance
- B: background error covariance

## Normalized Innovation
```
d_norm = d / √(σ²_total)
```

## Chi-Square Test
For N observations:
```
χ² = Σᵢ (d_i / σ_i)²
```
Expected value: χ² ≈ N if errors well-calibrated

## Decision Rule
Reject if |d_norm| > threshold where threshold is typically:
- Conservative: 3σ
- Standard: 5σ
- Liberal: 7σ (for sparse data)

# Platform-Specific Thresholds

**Radiosonde**: 3.0σ (high quality)
**Aircraft**: 3.5σ (temperature bias issues)
**Surface**: 2.5σ (high data density)
**Satellite**: 4.0σ (representation error)
**GPS-RO**: 3.0σ (bending angle), 4.0σ (refractivity)
"""
module BackgroundCheck

using LinearAlgebra
using Statistics
using ...ObservationTypes
using ...GSICoreAnalysis: AbstractStateVector

# Export types and functions
export BackgroundCheckConfig, BackgroundCheckResult
export perform_background_check, compute_innovation_statistics
export estimate_background_error, check_chi_square

"""
    BackgroundCheckConfig{T<:AbstractFloat}

Configuration for background check quality control.

# Fields
- `rejection_threshold::T`: Threshold for rejection (normalized units)
- `warning_threshold::T`: Threshold for warning (normalized units)
- `adaptive_threshold::Bool`: Use adaptive thresholds based on statistics
- `error_inflation_factor::T`: Factor to inflate obs errors for failed checks
- `use_ensemble_spread::Bool`: Use ensemble spread in background error
- `platform_thresholds::Dict{String,T}`: Platform-specific thresholds
- `variable_thresholds::Dict{ObservationVariable,T}`: Variable-specific thresholds
- `min_background_error::T`: Minimum background error (% of obs error)
"""
struct BackgroundCheckConfig{T<:AbstractFloat}
    rejection_threshold::T
    warning_threshold::T
    adaptive_threshold::Bool
    error_inflation_factor::T
    use_ensemble_spread::Bool
    platform_thresholds::Dict{String,T}
    variable_thresholds::Dict{ObservationVariable,T}
    min_background_error::T

    function BackgroundCheckConfig{T}(;
        rejection_threshold = T(3.0),
        warning_threshold = T(2.0),
        adaptive_threshold = true,
        error_inflation_factor = T(2.0),
        use_ensemble_spread = false,
        min_background_error = T(0.5)
    ) where T

        # Platform-specific thresholds
        platform_thresh = Dict{String,T}(
            "radiosonde" => T(3.0),
            "aircraft" => T(3.5),
            "surface" => T(2.5),
            "satellite" => T(4.0),
            "gps_ro" => T(3.5),
            "radar" => T(4.0),
            "ship" => T(3.0),
            "buoy" => T(3.0),
        )

        # Variable-specific thresholds
        variable_thresh = Dict{ObservationVariable,T}(
            TEMPERATURE => T(3.0),
            SURFACE_PRESSURE => T(2.5),
            WIND_U => T(3.5),
            WIND_V => T(3.5),
            SPECIFIC_HUMIDITY => T(4.0),
            BRIGHTNESS_TEMPERATURE => T(4.0),
            GPS_REFRACTIVITY => T(4.0),
            GPS_BENDING_ANGLE => T(3.0),
            RADAR_REFLECTIVITY => T(5.0),
            RADAR_VELOCITY => T(4.0),
        )

        new{T}(
            rejection_threshold,
            warning_threshold,
            adaptive_threshold,
            error_inflation_factor,
            use_ensemble_spread,
            platform_thresh,
            variable_thresh,
            min_background_error
        )
    end
end

BackgroundCheckConfig() = BackgroundCheckConfig{Float64}()

"""
    BackgroundCheckResult{T<:AbstractFloat}

Result of background check for a single observation.

# Fields
- `passed::Bool`: Whether observation passed background check
- `failure_reason::String`: Reason for failure (if failed)
- `innovation::T`: Innovation (obs - background)
- `normalized_innovation::T`: Normalized innovation
- `background_error::T`: Estimated background error
- `total_error::T`: Total error (obs + background)
- `qc_flag::Int`: QC flag (0=good, 1=suspect, 2=bad)
- `inflation_applied::Bool`: Whether error inflation was applied
"""
struct BackgroundCheckResult{T<:AbstractFloat}
    passed::Bool
    failure_reason::String
    innovation::T
    normalized_innovation::T
    background_error::T
    total_error::T
    qc_flag::Int
    inflation_applied::Bool
end

"""
    perform_background_check(obs::AbstractObservation{T},
                            background_value::T,
                            config::BackgroundCheckConfig{T};
                            background_error::Union{T,Nothing}=nothing,
                            ensemble_spread::Union{T,Nothing}=nothing) where T

Perform background check on observation.

# Arguments
- `obs`: Observation to check
- `background_value`: Background (first guess) value in observation space
- `config`: Background check configuration
- `background_error`: Optional background error estimate
- `ensemble_spread`: Optional ensemble spread for error estimation

# Returns
- `BackgroundCheckResult{T}`: Detailed check results

# Examples
```julia
config = BackgroundCheckConfig{Float64}()
result = perform_background_check(obs, H_xb, config)
if result.passed
    println("Innovation: ", result.innovation)
    println("Normalized: ", result.normalized_innovation)
else
    println("Rejected: ", result.failure_reason)
end
```
"""
function perform_background_check(obs::AbstractObservation{T},
                                 background_value::T,
                                 config::BackgroundCheckConfig{T};
                                 background_error::Union{T,Nothing}=nothing,
                                 ensemble_spread::Union{T,Nothing}=nothing) where T

    # Compute innovation
    innovation = obs.observed_value - background_value

    # Estimate background error if not provided
    if isnothing(background_error)
        background_error = estimate_background_error(obs, config, ensemble_spread)
    end

    # Compute total error (observation + background)
    obs_error = obs.observation_error.total_error
    total_error = sqrt(obs_error^2 + background_error^2)

    # Compute normalized innovation
    normalized_innovation = innovation / total_error

    # Get threshold for this observation
    threshold = get_threshold(obs, config)

    # Apply adaptive adjustment if enabled
    if config.adaptive_threshold
        threshold = adjust_threshold_adaptive(obs, threshold, config)
    end

    # Determine QC outcome
    abs_norm_innov = abs(normalized_innovation)

    passed = abs_norm_innov <= threshold
    inflation_applied = false

    qc_flag = if abs_norm_innov <= config.warning_threshold
        0  # Good
    elseif abs_norm_innov <= threshold
        1  # Suspect - passed but marginal
    else
        2  # Bad - failed check
    end

    failure_reason = if passed
        ""
    elseif abs_norm_innov <= threshold * 1.5
        "background_departure_moderate"
    else
        "background_departure_large"
    end

    return BackgroundCheckResult{T}(
        passed,
        failure_reason,
        innovation,
        normalized_innovation,
        background_error,
        total_error,
        qc_flag,
        inflation_applied
    )
end

"""
    estimate_background_error(obs::AbstractObservation{T},
                             config::BackgroundCheckConfig{T},
                             ensemble_spread::Union{T,Nothing}=nothing) where T

Estimate background error for observation location.

If ensemble spread is available and config enables it, use ensemble spread.
Otherwise, use observation error as proxy (typically 0.5-1.0 times obs error).
"""
function estimate_background_error(obs::AbstractObservation{T},
                                  config::BackgroundCheckConfig{T},
                                  ensemble_spread::Union{T,Nothing}=nothing) where T

    obs_error = obs.observation_error.total_error

    if config.use_ensemble_spread && !isnothing(ensemble_spread)
        # Use ensemble spread as background error estimate
        background_error = ensemble_spread
    else
        # Estimate based on observation type and conditions
        background_error = estimate_background_error_climatology(obs, obs_error)
    end

    # Ensure minimum background error
    min_bg_error = config.min_background_error * obs_error
    background_error = max(background_error, min_bg_error)

    return background_error
end

"""
    estimate_background_error_climatology(obs::AbstractObservation{T}, obs_error::T) where T

Estimate background error from climatology and observation characteristics.
"""
function estimate_background_error_climatology(obs::AbstractObservation{T}, obs_error::T) where T
    # Base estimate: typically 0.7-1.0 times observation error
    bg_error = T(0.8) * obs_error

    # Adjust based on observation type
    if obs isa RadiosondeObservation
        # Good vertical resolution - background error similar to obs error
        bg_error = T(0.9) * obs_error

        # Upper levels typically have larger background error
        if obs.pressure_level < 30000  # < 300 hPa
            bg_error *= T(1.2)
        end

    elseif obs isa SurfaceObservation
        # Dense network - background should match obs well
        bg_error = T(0.7) * obs_error

    elseif obs isa AircraftObservation
        # Aircraft temps have known biases - increase background error
        if obs.variable == TEMPERATURE
            bg_error = T(1.2) * obs_error
        end

    elseif obs isa RadianceObservation
        # Satellite radiances - representation error important
        bg_error = T(1.0) * obs_error

        # Cloud-affected scenes have larger error
        if obs.cloud_fraction > 0.3
            bg_error *= T(1.5)
        end

    elseif obs isa GPSROObservation
        # GPS-RO typically well-represented in model
        bg_error = T(0.8) * obs_error

        # Lower troposphere has more error
        if obs.location.altitude < 5000  # meters
            bg_error *= T(1.3)
        end

    elseif obs isa RadarObservation
        # Radar has large representation error
        bg_error = T(1.5) * obs_error
    end

    # Adjust for variable type
    if obs.variable == SPECIFIC_HUMIDITY || obs.variable == RELATIVE_HUMIDITY
        # Humidity has large background errors
        bg_error *= T(1.3)
    elseif obs.variable == PRECIPITATION_RATE
        # Precipitation is very difficult to forecast
        bg_error *= T(2.0)
    end

    return bg_error
end

"""
    get_threshold(obs::AbstractObservation{T}, config::BackgroundCheckConfig{T}) where T

Get background check threshold for observation.

Uses platform-specific and variable-specific thresholds if available.
"""
function get_threshold(obs::AbstractObservation{T}, config::BackgroundCheckConfig{T}) where T
    threshold = config.rejection_threshold

    # Check for platform-specific threshold
    platform_id = get_platform_type(obs)
    if haskey(config.platform_thresholds, platform_id)
        threshold = config.platform_thresholds[platform_id]
    end

    # Check for variable-specific threshold (takes precedence)
    if haskey(config.variable_thresholds, obs.variable)
        threshold = config.variable_thresholds[obs.variable]
    end

    return threshold
end

"""
    adjust_threshold_adaptive(obs::AbstractObservation{T}, threshold::T, config::BackgroundCheckConfig{T}) where T

Adaptively adjust threshold based on local conditions.
"""
function adjust_threshold_adaptive(obs::AbstractObservation{T}, threshold::T,
                                  config::BackgroundCheckConfig{T}) where T
    adjusted = threshold

    # Relax threshold in data-sparse regions
    # (In full implementation would check observation density)

    # Relax threshold for difficult-to-assimilate variables
    if obs.variable in [PRECIPITATION_RATE, CLOUD_FRACTION, RADAR_REFLECTIVITY]
        adjusted *= T(1.3)
    end

    # Tighten threshold for high-quality observations
    if obs isa RadiosondeObservation && obs.level_type == :mandatory
        adjusted *= T(0.9)
    end

    # Relax threshold near model boundaries
    if is_near_domain_boundary(obs)
        adjusted *= T(1.2)
    end

    return adjusted
end

"""
    compute_innovation_statistics(observations::Vector{<:AbstractObservation{T}},
                                 background_values::Vector{T},
                                 config::BackgroundCheckConfig{T}) where T

Compute innovation statistics for a set of observations.

Returns chi-square statistic, bias, RMS, and other diagnostics.
"""
function compute_innovation_statistics(observations::Vector{<:AbstractObservation{T}},
                                      background_values::Vector{T},
                                      config::BackgroundCheckConfig{T}) where T
    n_obs = length(observations)

    if n_obs == 0 || length(background_values) != n_obs
        return Dict{String,Any}()
    end

    # Compute innovations
    innovations = [obs.observed_value - background_values[i] for (i, obs) in enumerate(observations)]

    # Compute errors
    obs_errors = [obs.observation_error.total_error for obs in observations]
    bg_errors = [estimate_background_error(obs, config) for obs in observations]
    total_errors = sqrt.(obs_errors.^2 .+ bg_errors.^2)

    # Normalized innovations
    norm_innovations = innovations ./ total_errors

    # Chi-square statistic
    chi_square = sum(norm_innovations.^2)
    chi_square_per_obs = chi_square / n_obs

    # Statistics
    stats = Dict{String,Any}(
        "n_observations" => n_obs,
        "mean_innovation" => mean(innovations),
        "rms_innovation" => sqrt(mean(innovations.^2)),
        "std_innovation" => std(innovations),
        "mean_normalized_innovation" => mean(norm_innovations),
        "rms_normalized_innovation" => sqrt(mean(norm_innovations.^2)),
        "chi_square" => chi_square,
        "chi_square_per_obs" => chi_square_per_obs,
        "expected_chi_square" => Float64(n_obs),
        "chi_square_ratio" => chi_square / n_obs,
    )

    return stats
end

"""
    check_chi_square(statistics::Dict{String,Any}, tolerance::T=T(0.2)) where T

Check if chi-square statistic indicates well-calibrated errors.

Chi-square per observation should be close to 1.0 if errors are correct.
"""
function check_chi_square(statistics::Dict{String,Any}, tolerance::Float64=0.2)
    if !haskey(statistics, "chi_square_ratio")
        return false, "no_chi_square_statistic"
    end

    ratio = statistics["chi_square_ratio"]

    # Check if within tolerance
    if abs(ratio - 1.0) <= tolerance
        return true, "well_calibrated"
    elseif ratio < 1.0 - tolerance
        return false, "errors_overestimated"
    else
        return false, "errors_underestimated"
    end
end

# Helper functions

"""
    get_platform_type(obs::AbstractObservation) -> String

Get platform type identifier for observation.
"""
function get_platform_type(obs::SurfaceObservation)::String
    return "surface"
end

function get_platform_type(obs::RadiosondeObservation)::String
    return "radiosonde"
end

function get_platform_type(obs::AircraftObservation)::String
    return "aircraft"
end

function get_platform_type(obs::RadianceObservation)::String
    return "satellite"
end

function get_platform_type(obs::GPSROObservation)::String
    return "gps_ro"
end

function get_platform_type(obs::RadarObservation)::String
    return "radar"
end

function get_platform_type(obs::AbstractObservation)::String
    return "unknown"
end

"""
    is_near_domain_boundary(obs::AbstractObservation{T}) where T -> Bool

Check if observation is near model domain boundary.

(Simplified implementation - would need actual domain bounds)
"""
function is_near_domain_boundary(obs::AbstractObservation{T}) where T
    # Simplified: check for extreme latitudes
    if abs(obs.location.latitude) > 80.0
        return true
    end

    return false
end

end # module BackgroundCheck
